home *** CD-ROM | disk | FTP | other *** search
/ Plug-In Power Pack for Netscape Communicator / Plug-In Power Pack for Netscape Communicator.iso / plugins / dataviews / include / fdsevalfuns.h < prev    next >
C/C++ Source or Header  |  1997-05-08  |  29KB  |  956 lines
  1. #ifndef lint
  2. static char SccsId_FDSevalfuns[] = "@(#)FDSevalfuns.h    V1.32    3/16/95";
  3. #endif
  4. /*
  5. |    file name - FDSevalfuns.h
  6. |===================================================================
  7. |
  8. |              copyright (c) 1989,1992 V.I. Corporation
  9. |
  10. |    FDSevalfuns - Function definitions for FDSeval
  11. |
  12. |    Alan C Morse    31 may 89
  13. |    Alan C Morse    30 Dec 92     Added more functions: the
  14. |                    ANSI C functions we don't
  15. |                    yet support. atan2, ceil, exp,
  16. |                    fabs, floor, fmod, log, log10
  17. |                    Bit functions: bitand, bitor, bitxor
  18. |                    bitnot, bitget, bitset, bitclear;
  19. |                    Stat funs: runmax, runmedian,
  20. |                    runmean, runmin, and runsdev
  21. |    Alan C Morse    13 Jan 92    Added math funs pow and exp10
  22. |
  23. |===================================================================
  24. |
  25. |    Procedure description/function:
  26. |    Module containing function table and function definitions
  27. |    for FDSeval expression evaluator.
  28. |
  29. |    This is meant to be included in the eval.c module of FDSeval.
  30. |
  31. |    To add new functions:
  32. |        1) Write the function, type LOCAL double. If the function
  33. |        is a system function, you don't need a local function
  34. |        usually, but if you do, use the same name but with the
  35. |        first letter capitalized.
  36. |        2) Increment the NUMFCNS define
  37. |        3) Add it to FcnTable at the end of this file, giving
  38. |        3 fields: a) the name of the function as it would
  39. |        be typed in the expression; b) the name of the
  40. |        function in this module; and c) the number
  41. |        of arguments it expects. Putting things in
  42. |        alphabetical order and adding a comment field
  43. |        would be nice, too, but it is not required.
  44. |
  45. |===================================================================
  46. */
  47. #include "std.h"
  48. #include "vi/VUmiscfuns.h"
  49.  
  50. /*--------------------------------------------------------------*/
  51. /*             Functions for use in expressions            */
  52. /*--------------------------------------------------------------*/
  53. /* Simple functions; front ends for system function */
  54.  
  55. LOCAL double Abs( a )        double a;    { return DV_VIABS(a); }
  56.  
  57. /* -------------------------------------------------------
  58.     Exp10
  59.     exp10 is derived from the exp using the rule that
  60.         (a**b)**c : (a raised to b) raised to c
  61.         is the same as a**(b*c) : a raised to b*c
  62.     Hence:
  63.             10**a = (e ** log(10) ) ** a = e ** ( log10 * a )
  64. */
  65. LOCAL double NaturalLogOf10 = 0;
  66.  
  67. LOCAL double Exp10( a )
  68.   double a; 
  69.   {
  70.   if( NaturalLogOf10 == 0 )
  71.     NaturalLogOf10 = log(10.0);
  72.   return exp( a * NaturalLogOf10 );
  73.   }
  74.  
  75. LOCAL double Fmod( a, b )
  76.   double a,b; 
  77.   {
  78.   if( b != 0 )
  79.     return fmod(a,b);
  80.   else
  81.     return 0;
  82.   }
  83.  
  84. LOCAL double Log( a )
  85.   double a; 
  86.   {
  87.   if( a > 0 )
  88.     return log(a);
  89.   else
  90.     return 0;
  91.   }
  92.  
  93. LOCAL double Log10( a )
  94.   double a; 
  95.   {
  96.   if( a > 0 )
  97.     return log10(a);
  98.   else
  99.     return 0;
  100.   }
  101.  
  102. LOCAL double Min( a, b )    double a,b;    { return S_MIN(a,b); }
  103.  
  104. LOCAL double Max( a, b )    double a,b;    { return S_MAX(a,b); }
  105.  
  106. LOCAL double Pow( a, b )
  107.   double a,b; 
  108.   {
  109.   if( a > 0 )
  110.     return pow(a, b);
  111.   else
  112.     return 0;
  113.   }
  114.  
  115. /*==============================================================
  116. |                          BIT operations
  117. ================================================================*/
  118.  
  119. typedef unsigned long INT32BITS;
  120.  
  121. /*--------------------------------------------------------------
  122. |    bitand
  123. |    Returns bitwise AND of int32 versions of a and b.
  124. */
  125. LOCAL double bitand( a, b )
  126.   double a,b;
  127.   {
  128.   return (double)( (INT32BITS)a & (INT32BITS)b );
  129.   }
  130.  
  131. /*--------------------------------------------------------------
  132. |    bitclear
  133. |    CLEARs b-th bit of int32 version of a; 0 is least significant bit.
  134. |    Returns a after the operation.
  135. */
  136. LOCAL double bitclear( a, b )
  137.   double a,b;
  138.   {
  139.   INT32BITS ia, ib;
  140.  
  141.   ib = (INT32BITS)b;
  142.   ia = (INT32BITS)a;
  143.   /* Next line causes warning "unsigned value >= 0 is always 1".
  144.    * But keep it anyway in case definition of INT32BITS ever
  145.    * changes to be signed. 
  146.    */
  147.   if( ib >= 0 && ib <= 31 )    
  148.     ia = ~(1<<ib) & ia;
  149.   return (double)ia;
  150.   }
  151.  
  152. /*--------------------------------------------------------------
  153. |    bitget
  154. |    Returns the  b-th bit of int32 version of a; 0 is least significant bit.
  155. */
  156. LOCAL double bitget( a, b )
  157.   double a,b;
  158.   {
  159.   INT32BITS ia, ib;
  160.  
  161.   ib = (INT32BITS)b;
  162.   ia = (INT32BITS)a;
  163.   /* Next line causes warning "unsigned value >= 0 is always 1".
  164.    * But keep it anyway in case definition of INT32BITS ever
  165.    * changes to be signed. 
  166.    */
  167.   if( ib >= 0 && ib <= 31 )
  168.     if( ia & (1<<ib) )
  169.       return (double)1;
  170.     else
  171.       return (double)0;
  172.   else
  173.     return (double)0;
  174.   }
  175.  
  176. /*--------------------------------------------------------------
  177. |    bitnot
  178. |    Returns bitwise NOT of int32 version of a.
  179. */
  180. LOCAL double bitnot( a )
  181.   double a;
  182.   {
  183.   return (double)( ~(INT32BITS)a );
  184.   }
  185.  
  186. /*--------------------------------------------------------------
  187. |    bitor
  188. |    Returns bitwise OR of int32 versions of a and b.
  189. */
  190. LOCAL double bitor( a, b )
  191.   double a,b;
  192.   {
  193.   return (double)( (INT32BITS)a | (INT32BITS)b );
  194.   }
  195.  
  196. /*--------------------------------------------------------------
  197. |    bitset
  198. |    SETs b-th bit of int32 version of a; 0 is least significant bit.
  199. |    Returns a after the operation.
  200. */
  201. LOCAL double bitset( a, b )
  202.   double a,b;
  203.   {
  204.   INT32BITS ia, ib;
  205.  
  206.   ib = (INT32BITS)b;
  207.   ia = (INT32BITS)a;
  208.   /* Next line causes warning "unsigned value >= 0 is always 1".
  209.    * But keep it anyway in case definition of INT32BITS ever
  210.    * changes to be signed. 
  211.    */
  212.   if( ib >= 0 && ib <= 31 )
  213.     ia = (1<<ib) | ia;
  214.   return (double)ia;
  215.   }
  216.  
  217. /*--------------------------------------------------------------
  218. |    bitxor
  219. |    Returns bitwise XOR of int32 versions of a and b.
  220. */
  221. LOCAL double bitxor( a, b )
  222.   double a,b;
  223.   {
  224.   return (double)( (INT32BITS)a ^ (INT32BITS)b );
  225.   }
  226.  
  227. /*--------------------------------------------------------------
  228. |
  229. |    mk_normrand
  230. |    Makes a bunch of normally distributed random numbers.
  231. |    It fills a caller-supplied table with the numbers.
  232. |    Based on the Box-Muller transformation on p.453 of
  233. |    NUMERICAL METHODS.  The mean of the distribution is zero,
  234. |    and the standard deviation is one. Called by gnoise.
  235. */
  236.  
  237. LOCAL DV_BOOL L_randseeded = NO;
  238.  
  239. #ifndef PI
  240. #define PI      3.14159265358979323846264338327950
  241. #endif
  242.  
  243. LOCAL VOID mk_normrand( table, size )
  244.   FAST double *table;
  245.   int size;
  246.   {
  247.   int i;
  248.   double radius, angle;
  249.   double r1,r2,n1,n2;
  250.  
  251.   VUsrand(4321);
  252.   L_randseeded = YES;
  253.  
  254.   for( i=size/2; i>0; i-- )
  255.     {
  256.     r1 = VUfrand();
  257.     r2 = VUfrand();
  258.     if( r1 > 0 && r1 <= 1 )
  259.       radius = sqrt( -2.0 * log( r1 ) );
  260.     else
  261.       /* r1 is not valid, kludge up an alternate value */
  262.       /* If this happens, the distribution is no longer gaussian */
  263.       radius = DV_VIABS(r1);
  264.  
  265.     angle = 2 * PI * r2;
  266.     n1 = radius * cos( angle );
  267.     n2 = radius * sin( angle );
  268.     /* table entry = mean + standard_deviation * n[12] */
  269.     *table++ = n1;
  270.     *table++ = n2;
  271.     }
  272.   }
  273.  
  274. /*---------------------------------------------------------------- */
  275. /*    gnoise: Gaussian noise generator, w/ mean of 0 and stdev of 1 */
  276.  
  277. #define RAN_TABLE_SIZE    1024 /* Must be a power of two */
  278. LOCAL double *RanNumTable = 0;
  279.  
  280. LOCAL double gnoise()
  281.   {
  282.   if( RanNumTable == NULL )
  283.     {
  284.     RanNumTable = (double*)S_ALLOC( RAN_TABLE_SIZE * sizeof(double) );
  285.     mk_normrand( RanNumTable, RAN_TABLE_SIZE ); /* Fill it in */
  286.     }
  287.   return RanNumTable[ VUrand() & (RAN_TABLE_SIZE-1) ];
  288.   }
  289.   
  290. /*----------------------------------------------------------------
  291. |   Hypot:  return sqrt( a*a + b*b );
  292. */
  293. LOCAL double Hypot( a, b )
  294.      double a, b;
  295. {
  296.    return sqrt( a*a + b*b );
  297. }
  298.   
  299. /*---------------------------------------------------------------- */
  300. /*    wnoise:   white noise in range [0,1]            */
  301.  
  302. LOCAL double wnoise()
  303.  
  304.   {   
  305.   if (!L_randseeded) 
  306.     {
  307.     VUsrand(4321);
  308.     L_randseeded = YES;
  309.     }
  310.   return VUfrand();
  311.   }
  312.  
  313. /*---------------------------------------------------------------- */
  314. /*    saw:   Sawtooth wave generator.... /|/|/|/|        */
  315. /*     with data in range 0 to 1, with specified period    */
  316.  
  317. LOCAL double saw( iteration, period )
  318.   double iteration, period;
  319.   {
  320.   double p;
  321.  
  322.   /* Bring iteration into the range [0,period] */
  323.   p = iteration - ((int)(iteration/period)) *  period;
  324.  
  325.   return p / period;
  326.   }
  327.  
  328. /*---------------------------------------------------------------- */
  329. /*    rsaw:   Reverse sawtooth wave generator.... |\|\|\|\    */
  330. /*     with data in range 0 to 1, with specified period    */
  331.  
  332. LOCAL double rsaw( iteration, period )
  333.   double iteration, period;
  334.   {
  335.   double p;
  336.  
  337.   /* Bring iteration into the range [0,period] */
  338.   p = iteration - ((int)(iteration/period)) *  period;
  339.   p = period - p;
  340.  
  341.   return p / period;
  342.   }
  343.  
  344. /*---------------------------------------------------------------- */
  345. /*    tri:  Triangle wave generator.... /\/\/\/\                   */
  346. /*     with data in range 0 to 1, with specified period    */
  347.  
  348. LOCAL double tri( iteration, period )
  349.   double iteration, period;
  350.   {
  351.   double p;
  352.  
  353.   /* Bring iteration into the range [0,period] */
  354.   p = iteration - ((int)(iteration/period)) *  period;
  355.  
  356.   if( p > period/2.0 )
  357.     p = period - p;
  358.  
  359.   return 2 * p / period;
  360.   }
  361.  
  362. /*---------------------------------------------------------------- */
  363. /*    time:  Unix time function.  Returns the number of seconds    */
  364. /*    since Jan 1, 1970, 00:00.00.    */
  365. LOCAL double Time()
  366.   {
  367.   return (double)time( 0 );
  368.   }
  369.  
  370. /*---------------------------------------------------------------- */
  371. /*    delay:  Data stream a delayed by b time units.    */
  372. LOCAL VOID FreeDelayParams( p )
  373.   RING_BUFFER  *p;
  374.   {
  375.   if( p )
  376.     DESTROY_RING_BUFFER( p );
  377.   }
  378.  
  379. LOCAL double delay( a, delay_amount )
  380.   double a,delay_amount;
  381.   {
  382.   double *d;
  383.   int i, t;
  384.   double result;
  385.   RING_BUFFER *ringbuf; /* Local variable for ring buffer, since some
  386.     | compilers choke on CREATE_RING_BUFFER( GlobalEnv.CurrentNode->params... */
  387.  
  388.   ringbuf = (RING_BUFFER*)GlobalEnv.CurrentNode->params;
  389.   if( ringbuf == NULL )
  390.     { /* First call, create the ring buffer */
  391.     t = S_MAX( 1, (int)(delay_amount + 0.5) );
  392.     CREATE_RING_BUFFER( ringbuf, t, (LONG)sizeof(double) );
  393.     GlobalEnv.CurrentNode->params = (ADDRESS)ringbuf;
  394.     GlobalEnv.CurrentNode->FreeParams = FreeDelayParams;
  395.     
  396.     /* Prime the ring buffer with zeros */
  397.     for( i=0; i<t; i++ )
  398.       {
  399.       d = (double*)RING_BUFFER_ENTRY( ringbuf, i );
  400.       *d = (double)0;
  401.       }
  402.     }
  403.   /* Get the oldest value */
  404.   d = (double*)RING_BUFFER_ENTRY( ringbuf, 0 );
  405.   result = *d;
  406.   *d = a; /* Replace oldest value with newest value. */
  407.   INCREMENT_RING_BUFFER( ringbuf, 1 ); 
  408.  
  409.   return result;
  410.   }
  411.  
  412.  
  413. /*================================================================*/
  414. /*                 STATISTICS functions                           */
  415. /*                 over a time window                             */
  416. /*                 Prefixed by "run"                              */
  417. /*================================================================*/
  418.  
  419. /*----------------------------------------------------------------*/
  420. /*    Data buffer for time-windowed (running) functions.          */
  421. /*    The "RunningValue" datum is interpreted differently       */
  422. /*    depending on the function being performed.                */
  423. typedef struct
  424.   {
  425.   int iteration, width;
  426.   double RunningValue; /* Used for storing incremental intermediate data */
  427.   double RunningValue2; /* Used for storing incremental intermediate data */
  428.   RING_BUFFER *values; /* Buffered data values. There are width of them. */
  429.   } RUNDATA_PARAMS;
  430.  
  431. /* Allocate a RUNDATA_PARAMS structure */
  432. LOCAL RUNDATA_PARAMS *AllocRunParams( width )
  433.   double width;
  434.   {
  435.   int iwidth;
  436.   FAST RUNDATA_PARAMS *p;
  437.   
  438.   iwidth = S_MAX( 0, (int)(width+0.5) );
  439.   p = (RUNDATA_PARAMS*)S_ALLOC( sizeof(RUNDATA_PARAMS) );
  440.   p->iteration = 0;
  441.   p->width = S_MAX( 0, (int)(width+0.5) );
  442.   p->RunningValue = 0;
  443.   p->RunningValue2 = 0;
  444.   if( p->width )
  445.     { CREATE_RING_BUFFER( p->values, p->width, (LONG)sizeof(double) ); }
  446.   else
  447.     p->values = (RING_BUFFER*)0;
  448.   return p;
  449.   }
  450.  
  451. /* Free a RUNDATA_PARAMS structure */
  452. LOCAL VOID FreeRunParams( p )
  453.   RUNDATA_PARAMS *p;
  454.   {
  455.   if( p )
  456.     {
  457.     if( p->values ) DESTROY_RING_BUFFER( p->values );
  458.     S_FREE( (ADDRESS) p );
  459.     }
  460.   }
  461.  
  462. /*---------------------------------------------------------------- */
  463. /*    runmax:  The running max of a, over an interval of width points */
  464. /*    width should be a constant, since it is only examined the */
  465. /*    first time through.  If width is zero, the running max */
  466. /*    is over an infinitely wide period */
  467. LOCAL double runmax( a, width )
  468.   double a,width;
  469.   {
  470.   RUNDATA_PARAMS *p;
  471.   double *d;
  472.   DV_BOOL RecalcMax;
  473.   int i;
  474.  
  475.   if( GlobalEnv.CurrentNode->params == NULL )
  476.     {
  477.     GlobalEnv.CurrentNode->params = (ADDRESS)AllocRunParams( width );
  478.     GlobalEnv.CurrentNode->FreeParams = FreeRunParams;
  479.     p = (RUNDATA_PARAMS*)GlobalEnv.CurrentNode->params;
  480.     p->RunningValue = a; /* Initialize the maximum */
  481.     }
  482.   else
  483.     {
  484.     p = (RUNDATA_PARAMS*)GlobalEnv.CurrentNode->params;
  485.     if( a > p->RunningValue )
  486.       p->RunningValue = a;
  487.     }
  488.  
  489.   p->iteration++;
  490.   if( p->width )
  491.     {
  492.     d = (double*)RING_BUFFER_ENTRY( p->values, 0 ); /* Get oldest value */
  493.  
  494.     /* If we are about to delete the current max we will need to recalc. */
  495.     if( *d >= p->RunningValue )
  496.       RecalcMax = YES;
  497.     else
  498.       RecalcMax = NO;
  499.     *d = a; /* Replace oldest value with newest value. */
  500.     INCREMENT_RING_BUFFER( p->values, 1 ); /* Move past new value to next oldest */
  501.  
  502.     if( RecalcMax )
  503.       {
  504.       d = (double*)p->values->buffer;
  505.       p->RunningValue = *d++;
  506.       for( i = p->width-1; i>0; i--, d++ )
  507.     if( *d > p->RunningValue )
  508.       p->RunningValue = *d;
  509.       }
  510.     }
  511.  
  512.   return p->RunningValue;
  513.   }
  514.  
  515. /*---------------------------------------------------------------- */
  516. /*    runmin:  The running min of a, over an interval of width points */
  517. /*    width should be a constant, since it is only examined the */
  518. /*    first time through.  If width is zero, the running min */
  519. /*    is over an infinitely wide period */
  520. LOCAL double runmin( a, width )
  521.   double a,width;
  522. {
  523.   RUNDATA_PARAMS *p;
  524.   double *d;
  525.   DV_BOOL RecalcMin;
  526.   int i;
  527.  
  528.   if( GlobalEnv.CurrentNode->params == NULL )
  529.     {
  530.     GlobalEnv.CurrentNode->params = (ADDRESS)AllocRunParams( width );
  531.     GlobalEnv.CurrentNode->FreeParams = FreeRunParams;
  532.     p = (RUNDATA_PARAMS*)GlobalEnv.CurrentNode->params;
  533.     p->RunningValue = a; /* Initialize the minimum */
  534.     }
  535.   else
  536.     {
  537.     p = (RUNDATA_PARAMS*)GlobalEnv.CurrentNode->params;
  538.     if( a < p->RunningValue )
  539.       p->RunningValue = a;
  540.     }
  541.  
  542.   p->iteration++;
  543.   if( p->width )
  544.   {
  545.     d = (double*)RING_BUFFER_ENTRY( p->values, 0 ); /* Get oldest value */
  546.  
  547.     /* If we are about to delete the current min we will need to recalc. */
  548.     if( p->iteration > 1 && *d <= p->RunningValue )
  549.       RecalcMin = YES;
  550.     else
  551.       RecalcMin = NO;
  552.     *d = a; /* Replace oldest value with newest value. */
  553.     INCREMENT_RING_BUFFER( p->values, 1 ); /* Move past new value to next oldest */
  554.  
  555.     if( RecalcMin )
  556.     {
  557.       d = (double*)p->values->buffer;
  558.       p->RunningValue = *d++;
  559.       for( i = S_MIN(p->width-1,p->iteration-1); i>0; i--, d++ )
  560.                 if( *d < p->RunningValue )
  561.                     p->RunningValue = *d;
  562.         }
  563.     }
  564.  
  565.   return p->RunningValue;
  566. }
  567.  
  568. /*---------------------------------------------------------------- */
  569. /*    runmean:  The running mean of a, over an interval of width points */
  570. /*    width should be a constant, since it is only examined the */
  571. /*    first time through.  If width is zero, the running mean */
  572. /*    is over an infinitely wide period */
  573. LOCAL double runmean( a, width )
  574.   double a,width;
  575. {
  576.   RUNDATA_PARAMS *p;
  577.   double *d;
  578.   double result;
  579.  
  580.   if( GlobalEnv.CurrentNode->params == NULL )
  581.     {
  582.     GlobalEnv.CurrentNode->params = (ADDRESS)AllocRunParams( width );
  583.     GlobalEnv.CurrentNode->FreeParams = FreeRunParams;
  584.     }
  585.   p = (RUNDATA_PARAMS*)GlobalEnv.CurrentNode->params;
  586.   p->iteration++;
  587.   p->RunningValue += a;
  588.   if( p->width )
  589.     {
  590.     d = (double*)RING_BUFFER_ENTRY( p->values, 0 ); /* Get oldest value */
  591.     if( p->iteration > p->width )
  592.         { p->RunningValue -= *d; result = p->RunningValue / (double)p->width; }
  593.     else
  594.       result = p->RunningValue / (double)p->iteration;
  595.     *d = a; /* Replace oldest value with newest value. */
  596.     INCREMENT_RING_BUFFER( p->values, 1 ); /* Move past new value to next oldest */
  597.     }
  598.   else
  599.     result = p->RunningValue / (double)p->iteration;
  600.  
  601.   return result;
  602. }
  603.  
  604. /*---------------------------------------------------------------- */
  605. /*    runmedian:  The running median of a, over an interval of width points */
  606. /*    width should be a constant, since it is only examined the */
  607. /*    first time through.  Unlike the mean, if width is zero,   */
  608. /*    the running median can't be calculated over an infinite   */
  609. /*      because the median calculation needs all the data to be saved,  */
  610. /*    which would be too expensive.  */
  611.  
  612. /* Iterative algorithm for finding the median of a data set without sorting */
  613. /* If you pass in a reasonable first guess than the algorithm converges faster*/
  614. /* This algorighm comes from Numerical Recipes In C, page 477-479. */
  615.  
  616. /* The median function needs DBL_MAX and DBL_MIN which are supposed to */
  617. /* be in float.h. The sun doesn't have this include file in SunOS, but */
  618. /* it does have values.h, which has MAXDOUBLE and MINDOUBLE */
  619. #ifndef DBL_MAX
  620. #  ifdef SUNOS
  621. #    include <values.h>
  622. #    define DBL_MAX MAXDOUBLE
  623. #    define DBL_MIN MINDOUBLE
  624. #  else
  625. #    ifdef MASSCOMP
  626. #      include <values.h>
  627. #      define DBL_MAX MAXDOUBLE
  628. #      define DBL_MIN MINDOUBLE
  629. #    else
  630. #    include <float.h>
  631. #  endif
  632. #endif
  633. #endif
  634. /* Amplification factors used in median calculation */
  635. #define MEDIAN_AFAC 1.5
  636. #define MEDIAN_AMP 1.5
  637.  
  638. LOCAL double median(x, n, FirstGuess)
  639.   double *x; /* Array of data values to find median */
  640.   int n; /* Size of array of data */
  641.   double *FirstGuess; /* First guess of median value;
  642.              if NULL make a guess in the routine. */
  643. {
  644.   int NumAbove, /* Number of points above the current guess */
  645.         NumBelow, /* Number of points below the current guess */
  646.         i;
  647.   double xmed; /* Median */
  648.   double Guess, /* Current Guess at the median */
  649.         NextGuess, /* Next Guess at the median (temporary variable) */
  650.         GuessUpperBound, /* Upper bound on the median. */
  651.         GuessLowerBound, /* Lower bound on the median. */
  652.         AboveGuess, /* Value above guess and closest to it in the data set. */
  653.         BelowGuess, /* Value below guess and closest to it in the data set. */
  654.         Error; /* Difference between guess and calculated median. */
  655.   FAST double *datum;
  656.   double temp, sumx, sum, eps;
  657.   
  658.   /* Make a guess at the median */
  659.   if( FirstGuess )
  660.     Guess = *FirstGuess;
  661.   else
  662.     Guess = 0.5 * (x[0] + x[n-1]);
  663.   eps = fabs(x[n-1] - x[0]);
  664.  
  665.   /* First guess at characteristinc spacing of the data points near the median */
  666.   GuessLowerBound = -DBL_MAX;
  667.   GuessUpperBound = DBL_MAX;
  668.   
  669.   FOREVER
  670.     { /* Start a pass through the data */
  671.             sum = sumx = 0.0;
  672.             NumAbove = NumBelow = 0;
  673.             BelowGuess = -DBL_MAX;
  674.             AboveGuess = DBL_MAX;
  675.  
  676.             /* For each data point */
  677.             for( datum = x, i = n; i > 0; datum++, i-- )
  678.       {
  679.                 if (*datum != Guess)
  680.                 {
  681.                     if (*datum > Guess)
  682.                     {
  683.                         ++NumAbove;
  684.                         if (*datum < AboveGuess)
  685.                             AboveGuess = *datum;
  686.                     }
  687.                     else if (*datum < Guess)
  688.                     {
  689.                         ++NumBelow;
  690.                         if (*datum > BelowGuess)
  691.                             BelowGuess = *datum;
  692.                     }
  693.                     sum += temp = 1.0 / (eps + fabs(*datum - Guess));
  694.                     sumx += *datum * temp;
  695.                 }
  696.       }
  697.             Error = (sumx / sum) - Guess;
  698.             if (NumAbove - NumBelow >= 2)
  699.       {
  700.                 GuessLowerBound = Guess;
  701.                 NextGuess = Error < 0.0 ? AboveGuess : AboveGuess + Error * MEDIAN_AMP;
  702.                 if (NextGuess > GuessUpperBound)
  703.                     NextGuess = 0.5 * (Guess + GuessUpperBound);
  704.                 eps = MEDIAN_AFAC * fabs(NextGuess - Guess);
  705.                 Guess = NextGuess;
  706.       }
  707.             else if (NumBelow - NumAbove >= 2)
  708.       {
  709.                 GuessUpperBound = Guess;
  710.                 NextGuess = Error > 0.0 ? BelowGuess : BelowGuess + Error * MEDIAN_AMP;
  711.                 if (NextGuess < GuessLowerBound)
  712.                     NextGuess = 0.5 * (Guess + GuessLowerBound);
  713.                 eps = MEDIAN_AFAC * fabs(NextGuess - Guess);
  714.                 Guess = NextGuess;
  715.       }
  716.             else
  717.       { /* Got the median */
  718.                 if (n % 2 == 0)
  719.                 { /* For even n, median is always an average. */
  720.                     if( NumAbove == NumBelow )
  721.                         xmed = AboveGuess + BelowGuess;
  722.                     else if (NumAbove > NumBelow )
  723.                         xmed = Guess + AboveGuess;
  724.                     else
  725.                         xmed = BelowGuess + Guess;
  726.                     xmed *= 0.5;
  727.                 }
  728.                 else
  729.                 { /* For odd n, median is always one point. */
  730.                     if( NumAbove == NumBelow )
  731.                         xmed = Guess;
  732.                     else if( NumAbove > NumBelow )
  733.                         xmed = AboveGuess;
  734.                     else
  735.                         xmed = BelowGuess;
  736.                 }
  737.                 return xmed;
  738.       }
  739.     }
  740. }
  741. /*----------------------------------------------------------*/
  742. LOCAL double runmedian( a, width )
  743.   double a,width;
  744. {
  745.   RUNDATA_PARAMS *p;
  746.   double *d;
  747.  
  748.   if( GlobalEnv.CurrentNode->params == NULL )
  749.     {
  750.     if( width <= 0 ) width = 100;
  751.     GlobalEnv.CurrentNode->params = (ADDRESS)AllocRunParams( width );
  752.     GlobalEnv.CurrentNode->FreeParams = FreeRunParams;
  753.     }
  754.   p = (RUNDATA_PARAMS*)GlobalEnv.CurrentNode->params;
  755.   p->iteration++;
  756.  
  757.   if( p->iteration <= p->width )
  758.     {
  759.     /* Add the new value to the data buffer at next empty slot. */
  760.     d = (double*)RING_BUFFER_ENTRY( p->values, p->iteration-1 );
  761.     *d = a;
  762.  
  763.     /* Recalculate the median */
  764.     p->RunningValue = median( (double*)p->values->buffer,
  765.                                                             p->iteration, &p->RunningValue );
  766.     }
  767.   else
  768.     {
  769.     /* Replace the oldest value with the current value */
  770.     d = (double*)RING_BUFFER_ENTRY( p->values, 0 ); /* Get oldest value */
  771.     *d = a;
  772.  
  773.     /* Recalculate the median */
  774.     p->RunningValue = median( (double*)p->values->buffer,
  775.                                                             p->width, &p->RunningValue );
  776.  
  777.     /* Increment the ring buffer to move past new value to next oldest */
  778.     INCREMENT_RING_BUFFER( p->values, 1 );
  779.  
  780.     }
  781.  
  782.   return p->RunningValue;
  783. }
  784.  
  785. /*---------------------------------------------------------------- */
  786. /*    runsdev:  The running standard deviation of a, over an interval */
  787. /*      of width points. Width should be a constant, since it is only */
  788. /*    examined the first time through.  If width is zero, the running */
  789. /*    standard deviation is over an infinitely wide period */
  790.  
  791. /* The equation that is used is
  792.    sqrt( ( sum of squares(a) - square of sum(a) / n )  / ( n - 1 ) )
  793.  
  794.    This algorithm is not computationally optimal (according to,
  795.    Numerical Methods in C, pg 475), but it is suitable for incremental
  796.    (one-pass) calculation.
  797. */
  798. LOCAL double runsdev( a, width )
  799.   double a,width;
  800. {
  801.   RUNDATA_PARAMS *p;
  802.   double *d;
  803.   double result, num;
  804.  
  805.   if( GlobalEnv.CurrentNode->params == NULL )
  806.     {
  807.     GlobalEnv.CurrentNode->params = (ADDRESS)AllocRunParams( width );
  808.     GlobalEnv.CurrentNode->FreeParams = FreeRunParams;
  809.     }
  810.   p = (RUNDATA_PARAMS*)GlobalEnv.CurrentNode->params;
  811.   p->iteration++;
  812.   p->RunningValue += a; /* Running sum */
  813.   p->RunningValue2 += a * a; /* Running sum of squares. */
  814.   if( p->width )
  815.     {
  816.     d = (double*)RING_BUFFER_ENTRY( p->values, 0 ); /* Get oldest value */
  817.     if( p->iteration > p->width )
  818.         {
  819.       p->RunningValue -= *d;
  820.       p->RunningValue2 -= ((*d) * (*d));
  821.       num = p->width;
  822.         }
  823.     else
  824.       num = p->iteration;
  825.     *d = a; /* Replace oldest value with newest value. */
  826.     INCREMENT_RING_BUFFER( p->values, 1 ); /* Move past new value to next oldest */
  827.     }
  828.   else
  829.     num = p->iteration;
  830.  
  831.   if( num > 1 )
  832.     result = sqrt( ( p->RunningValue2 -
  833.                                          p->RunningValue * p->RunningValue / num ) / ( num - 1 ) );
  834.   else
  835.     result = 0;
  836.    
  837.   return result;
  838. }
  839.  
  840. /*---------------------------------------------------------------- */
  841. /*    delta:  The change per unit time in a */
  842. /*          averaged over the last time_span time units */
  843. /*     NOTE: Use mean's data structure for storing parameters */
  844. #define MAX_TIME_SPAN 2048
  845.  
  846. typedef struct
  847. {
  848.   int iteration, width;
  849.   double FirstValue;
  850.   RING_BUFFER *values;
  851. } DELTA_PARAMS;
  852.  
  853. LOCAL VOID FreeDeltaParams( p )
  854.   DELTA_PARAMS *p;
  855. {
  856.   if( p )
  857.     {
  858.     if( p->values ) DESTROY_RING_BUFFER( p->values );
  859.     S_FREE( (ADDRESS) p );
  860.     }
  861. }
  862.  
  863. LOCAL double delta( a, time_span )
  864.   double a, time_span;
  865. {
  866.   DELTA_PARAMS *p;
  867.   double *d;
  868.   double result;
  869.  
  870.   p = (DELTA_PARAMS*)GlobalEnv.CurrentNode->params;
  871.   if( GlobalEnv.CurrentNode->params == NULL )
  872.     {
  873.     GlobalEnv.CurrentNode->params = S_ALLOC( sizeof(DELTA_PARAMS) );
  874.     p = (DELTA_PARAMS*)GlobalEnv.CurrentNode->params;
  875.     p->iteration = 0;
  876.     p->width = S_MIN( MAX_TIME_SPAN, S_MAX( 1, (int)(time_span+0.5) ) );
  877.     p->FirstValue = a;
  878.     if( p->width )
  879.         { CREATE_RING_BUFFER( p->values, p->width, (LONG)sizeof(double) ); }
  880.  
  881.     GlobalEnv.CurrentNode->FreeParams = FreeDeltaParams;
  882.     }
  883.   p->iteration++;
  884.   d = (double*)RING_BUFFER_ENTRY( p->values, 0 ); /* Get oldest value */
  885.  
  886.   if( p->iteration > p->width )
  887.     result = (a - *d) / (double)p->width;
  888.   else
  889.     /* Haven't run enough iterations yet */
  890.     /* Assume oldest value is is the same as the first value */
  891.     result = (a - p->FirstValue) / (double)p->iteration;
  892.  
  893.   *d = a; /* Replace oldest value with newest value. */
  894.  
  895.   /* Move past new value to next oldest */
  896.   INCREMENT_RING_BUFFER( p->values, 1 );
  897.  
  898.   return result;
  899. }
  900.  
  901. /*--------------------------------------------------------------*/
  902. /*             FUNCTION TABLE                */
  903. /*--------------------------------------------------------------*/
  904.  
  905. /* The typedef for FUN_DESC is in FDSeval.c */
  906.  
  907. #define NUMFCNS 45
  908. LOCAL FUN_DESC FcnTable[NUMFCNS] =
  909. { /* Note for comments, arguments are (a,b,c...) */
  910.   "abs",    Abs,    1, /* absolute value of a */
  911.   "acos",    acos,    1, /* arc cosine of a */
  912.   "asin",    asin,    1, /* arc sine of a */
  913.   "atan",    atan,    1, /* arc tangent of a */
  914.   "atan2",    atan2,    2, /* arc tangent of a/b */
  915.   "bitand",    bitand, 2, /* bitwise AND of int32 versions of a and b */
  916.   "bitclear",    bitclear,2, /* CLEARs b-th bit of int32 version of a; 0 is lsb */
  917.   "bitget",    bitget,    2, /* gets b-th bit of int32 version of a; 0 is lsb */
  918.   "bitnot",    bitnot,    1, /* bitwise NOT of int32 version of a */
  919.   "bitor",    bitor,     2, /* bitwise OR of int32 versions of a and b */
  920.   "bitset",    bitset,    2, /* SETs b-th bit of int32 version of a; 0 is lsb */
  921.   "bitxor",    bitxor,    2, /* bitwise XOR of int32 versions of a and b */
  922.   "ceil",    ceil,    1, /* smallest integer not less than a */
  923.   "cos",    cos,    1, /* cosine of a */
  924.   "cosh",    cosh,    1, /* hyperbolic cosine of a */
  925.   "delay",    delay,    2, /* a delayed by b iterations */
  926.   "delta",    delta,    2, /* change in a averaged over b iterations */
  927.   "exp",    exp,    1, /* exponential of a; i.e. e raised to the a power */
  928.   "exp10",    Exp10,    1, /* base 10 exponential of a; i.e. 10 raised to the a power */
  929.   "fabs",    fabs,    1, /* absolute value of a (same as DV_VIABS) */
  930.   "floor",    floor,    1, /* largest integer value not greater than a */
  931.   "fmod",    Fmod,    2, /* remainder of a/b */
  932.   "gnoise",    gnoise,    0, /* gaussian noise, with mean 0, std dev 1 */
  933.   "hypot",    Hypot,    2, /* sqrt( a*a + b*b ) */
  934.   "log",    Log,    1, /* the log base e of a */
  935.   "log10",    Log10,    1, /* the log base 10 of a */
  936.   "max",    Max,    2, /* the maximum of a and b */
  937.   "mean",    runmean,2, /* running mean of a over a b-width window */
  938.   "min",    Min,    2, /* the minimum of a and b */
  939.   "pow",    Pow,    2, /* a raised to the floor(b) power; a > 0 */
  940.   "rsaw",    rsaw,    2, /* a=iteration, b=period of reverse saw */
  941.   "runmax",    runmax, 2, /* running max of a over a b-width window */
  942.   "runmean",    runmean,2, /* running mean of a over a b-width window */
  943.   "runmedian",    runmedian,2, /* running median of a over a b-width window (b>0) */
  944.   "runmin",    runmin, 2, /* running min of a over a b-width window */
  945.   "runsdev",    runsdev,2, /* running standard deviation of a over a b-width window */
  946.   "saw",    saw,    2, /* a=iteration, b=period of sawtooth*/
  947.   "sin",    sin,    1, /* sine of a */
  948.   "sinh",    sinh,    1, /* hyperbolic sin of a */
  949.   "sqrt",    sqrt,    1, /* square root of a */
  950.   "tan",    tan,     1, /* tangent of a */
  951.   "tanh",    tanh,     1, /* hyperbolic tangent of a */
  952.   "time",    Time,     0, /* number of seconds since 1 jan 1970 */
  953.   "tri",    tri,     2, /* a=iteration, b=period of triangle wave */
  954.   "wnoise",    wnoise,    0, /* white noise in the range [0,1] */
  955. };
  956.